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O ■ Abstract 
(N 

, Axisymmetric numerical simulations of rotating stellar core collapse to a 

neutron star are performed in the framework of full general relativity. The 
>^ , so-called Cartoon method, in which the Einstein field equations are solved in 

o . 

tJ" , the Cartesian coordinates and the axisymmetric condition is imposed around 

o 

O 

i cylindrical coordinates (on the y = plane in the Cartesian coordinates) 

using a high-resolution shock-capturing scheme with the maximum grid size 
(2500,2500). A parametric equation of state is adopted to model collapsing 



the y = plane, is adopted. The hydrodynamic equations are solved in the 



stellar cores and neutron stars following Dimmelmeier et al. It is found that 
the evolution of central density during the collapse, bounce, and formation 
of protoneutron stars agree well with those in the work of Dimmelmeier et 
al. in which an approximate general relativistic formulation is adopted. This 
indicates that such approximation is appropriate for following axisymmetric 
stellar core collapses and subsequent formation of protoneutron stars. Gravi- 
tational waves are computed using a quadrupole formula. It is found that the 
waveforms are qualitatively in good agreement with those by Dimmelmeier 
et al. However, quantitatively, two waveforms do not agree well. Possible 
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reasons for the disagreement are discussed. 
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I. INTRODUCTION 



Rotating stellar core collapse is among the most promising sources of gravitational waves. 
To date, there has been no systematic work for computation of the collapse to a neutron 
star and of emitted gravitational waves in full general relativity (but see [1]). Gravitational 
waves associated with the formation of rotating neutron stars have been widely computed 
in the Newtonian gravity [2-9] or in an approximate general relativistic gravity [10] using 
the so-called conformal flatness approximation (or Isenberg- Wilson-Mathews approximation 
[11]). As indicated in [10], the general relativistic effects modify the dynamics of the collapse 
and corresponding gravitational waveforms significantly. This implies that simulation in full 
general relativity is the best approach for accurate computation of gravitational waves. 

During the stellar core collapse to a neutron star, the characteristic radius changes from 
initial stellar core radius ~ 2000 km to neutron star radius ~ 10 km. Adopting a uniform 
and fixed grid with a grid spacing as ~ 1 km, the required grid number for the simulation is 
more than 2000 for one direction. With current computational resources, it is very difficult 
to take such a huge number of grid points in three-dimensional simulations. If the progenitor 
of the neutron star is not very rapidly rotating, nonaxisymmetric instabilities will not set 
in and the collapse will proceed in an axisymmetric manner. By restricting our attention 
to axisymmetric spacetimes, the grid resolution can be improved significantly for a given 
computational resource. Thus, as a first step, it is better to perform axisymmetric simula- 
tions than to do nonaxisymmetric ones for a well-resolved and convergent computation of 
the collapse, the bounce, and corresponding gravitational waveforms, focusing only on the 
moderately rapid rotation case. 

In this paper, we study gravitational waves from rotating stellar core collapses to a 
neutron star assuming the axial symmetry. Dynamics of the collapse is followed by fully 
general relativistic simulation. Gravitational waves are approximately computed using a 
quadrupole formula adopted and tested in [12]. Necessity of adopting quadrupole formulas 
arises from the fact that the amphtude of gravitational waves is too small (< 10""^ in the 
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local wave zone) to accurately extract the waveforms from raw data sets of metric. Although 
exact gravitational waveforms cannot be computed, the quadrupole formula is a useful tool 
for approximate computation of gravitational waves associated with matter motion such as 
oscillations of neutron stars as indicated in [12]. 

Recently, gravitational waves from axisymmetric rotating stellar core collapses have been 
extensively computed in a relativistic manner by Dimmelmeier et al. [10]. As mentioned 
above, they determine the gravitational fields adopting an approximate formulation of the 
Einstein equation. The approximation is likely to be applicable to a moderately relativistic 
and stationary spacetime such as that for a rapidly rotating neutron star [13]. However, no 
one has clarified whether this is the case for dynamical spacetimes. To confirm that their 
treatment is indeed appropriate, it is necessary to compare their solutions with fully general 
relativistic ones for a calibration. One of the purposes in this paper is to examine whether 
the numerical solution for the stellar core collapse computed in [10] is a well approximated 
one for a fully general relativistic solution. 

In [10], gravitational waveforms were computed in terms of a quadrupole formula. In gen- 
eral relativity, there is no unique definition of the quadrupole moment, nor is the quadrupole 
formula, for axisymmetric dynamical spacetimes. Accuracy of gravitational waveforms de- 
pends on the choice of the quadrupole formula and the gauge conditions. Thus, to know 
how accurately the approximate gravitational waveforms can be computed by the chosen 
quadrupole formula, a calibration is required by comparing the resulting waveforms with 
those computed from the metric as we did in [12]. Unfortunately, such calibration is not 
possible in [10] since they adopted an approximate general relativistic formulation for the 
gravitational field in which the metric does not carry any information of gravitational waves. 
Consequently, it is not clear whether the quadrupole formula they adopted can actually yield 
accurate approximate gravitational waveforms and how large the magnitude of the error is. 
On the other hand, in the previous paper [12], we did such a calibration for a quadrupole 
formula which is different from that in [10], and showed it possible to compute gravitational 
waves from oscillating and rapidly rotating neutron stars of high values of compactness 
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fairly accurately, besides possible systematic errors for the amplitude due to neglecting post- 
Newtonian corrections. Computing gravitational waveforms by the calibrated quadrupole 
formula and comparing the results with the previous ones, we estimate how accurate the 
waveforms computed in [10]. 

This paper is organized as follows. In Sec. II, our numerical implementations for general 
relativistic simulation in axial symmetry are briefly reviewed. In Sec. Ill, the initial condition 
and computational setting are described. Sec. IV presents the numerical results. Sec. V 
is devoted to a summary. Throughout this paper, we adopt the geometrical units in which 
G = c = 1 where G and c are the gravitational constant and the speed of light, respectively. 



II. NUMERICAL IMPLEMENTATION 

A. Summary of formulation 

We perform fully general relativistic simulations for rotating stellar core collapse in axial 
symmetry using the same formulation as in [14], to which the reader may refer for details 
and basic equations. The fundamental variables for the hydrodynamics are 

p: rest mass density, 
e: specific internal energy, 
P: pressure, 
u^: four velocity, 

where subscripts i,j,k,--- denote x,y and z, and /i the spacetime components. As the 
variables to be evolved in the numerical simulations, we define a weighted density p*(= 
pau^e^'^) and a weighted four- velocity Ui[= (1 -|- £ -|- P/p)ui\. From these variables, the 
total baryon rest-mass and angular momentum of system, which are conserved quantities in 
axisymmetric spacetimes, can be defined as 
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r* = y d^xp^, (2) 
J = j d^xpS^. (3) 



General relativistic hydrodynamic equations are solved using a so-called high-resolution 
shock-capturing scheme [15,14] on the y = plane with the cylindrical coordinates (x, z) (in 
the Cartesian coordinates with y = 0). 

The fundamental variables for geometry are 

a: lapse function, 
(3^: shift vector, 

7ij: metric in three dimensional spatial hypersurface, 
7= ei2^ = det(7,,), 

Kij-: extrinsic curvature. (4) 

We evolve 7ij, 0, 74jj = e~^'^{Kij — ^ijK^'^), and trace of the extrinsic curvature Kf} together 
with three auxiliary functions = S^'^dj'^ik with an unconstrained free evolution code as in 
[16,18-20,17,14]. 

The Einstein equations are solved in the Cartesian coordinates. To impose axisymmetric 
boundary conditions, the Cartoon method is used [22]: Assuming the reflection symmetry 
with respect to the equatorial plane, simulations are performed using a fixed uniform grid 
with the grid size N x3 x N in {x,y, z) which covers a computational domain as < x < L, 
< z < L, and —A < y < A. Here, and L are constants and A = L/N. In the Cartoon 
method, the axisymmetric boundary conditions are imposed at y — ±A. 

As the time slice, we impose an "approximate" maximal slicing condition in which K^J^ 
is required [16]. As the spatial gauge, we adopt a dynamical gauge condition [21] in which 
the equation for the shift vector is written as 

dtP'^^'^iFi + AtdtFi), (5) 



where At denotes a timestep in numerical computation. 

During the numerical simulations, violations of the Hamiltonian constraint and conser- 
vation of mass and angular momentum are monitored as code checks. Numerical results for 
several test calculations, including stability and collapse of nonrotating and rotating neutron 
stars, have been described in [14]. 

An outgoing wave boundary condition for Fj, hij{= jij — 5ij), and Aij is imposed at the 
outer boundaries of the computational domain. The condition adopted is the same as that 
described in [20]. K),'' is set to be zero at the outer boundaries. 



A parametric equation of state is adopted following Miiller and his collaborators [6,10]. 
In this equation of state, one assumes that the pressure consists of the sum of polytropic 
and thermal parts as 



The polytropic part is given by Pp = Kp{p)p^^P^ where Kp and F are not constants but 
functions of p. In this paper, we follow [10] for the choice of Kp{p) and r(p): For the density 
smaller than the nuclear density which is defined as pnuc = 2 X 10^^ g/cm^, F = Fi(=const) 
is set to be ^ 4/3, and for p > pnuc, T = F2(= const) > 2. Thus, 



[ K2P \ P> Pnuc, 

where Ki and K2 are constants. Since Pp should be continuous, the relation, K2 = Kip^^~^^, 
is required. Following [6,10], the value of Ki is fixed to 5 x 10^^ cgs. With this choice, the 
polytropic part of the equation of state for p < p^^^, in which the degenerate pressure of 
electrons is dominant, is approximated well. Since the specific internal energy should be 
continuous at p = pnuo the polytropic specific internal energy ep is defined as 



B. Equations of state 



P = Pp + Pth. 



(6) 
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+ 



(r^ - TiWn^- 
(r,-i)(r2-i) 



, p > p. 



p < Pi 



nuc- 



(8) 



With this setting, a reahstic equation of state for high-density, cold nuclear matter is mim- 
icked. 

The thermal part of the pressure Pth plays an important role in the case that shocks are 
generated. Pth is related to the thermal energy density Sth = £ — £p as 



Following [10], the value of Fth, which determines the strength of shocks, is chosen as 1.5 
for most simulations in this paper. Extending the previous work [10], for a few models, we 
set Fth = 1.35 or 5/3 to investigate the effect of the shock heating at the bounce phase and 
resulting gravitational waveforms. 

Simulations are initiated in the following manner: First, equilibrium rotating stars with 
F = 4/3 polytrope are given. Then, to induce the collapse, we shghtly decrease the value of 
the adiabatic index from F = 4/3 to Fi < 4/3. The equilibrium states are computed with 
the polytropic equation of state as 



where following [10], Kq is set to be 5 x 10^^ cm^/s^/gr^/^, with which a soft equation of 
state governed by the electron degenerate pressure is approximated well [23]. Here, Kq and 



-Pth = (rth - i)p£th- 



(9) 



P = Kop'/\ 



(10) 



Ki are related by Ki = KqPq 



4/3-r 



where we set po = I g/cm^. 



C. Quadrupole formula 



In the present work, gravitational waveforms are computed using a quadrupole formula 
[12]. In quadrupole formulas, gravitational waves at null infinity are calculated from 



where and Pi'' — — riin^ (n* = x^r) denote a tracefree quadrupole moment and the 
projection tensor. Prom this expression, the +-mode of gravitational waves with Z = 2 in 
axisymmetric spacetimes is written as 

where lij denotes a quadrupole moment, /y its second time derivative, and tret a retarded 
time. 

In fully general relativistic and dynamical spacetimes, there is no unique definition for 
the quadrupole moment and nor is for 7^. Following a previous paper [12], we choose the 
simplest definition as 

!■■ = j p^x'x^d^x. (13) 

Then, using the continuity equation of the form 

dtp, + diip^) ^ 0, (14) 

the first time derivative can be written as 

iij ^ j p^{v'x^ + x'v^)d^x. (15) 

To compute 7^, the finite differencing of the numerical result for 7^ is carried out. 

As indicated in [12], it is possible to compute gravitational waves from oscillating and 
rapidly rotating neutron stars of high values of compactness fairly accurately with the present 
choice of 7jj , besides possible systematic errors for the amphtude of order M/ R or v'^/ (? where 
M, 7?, and v denote the gravitational mass, the equatorial circumferential radius, and the 
radial velocity of the collapsing star and/or formed neutron stars. In stellar core collapses, 
l(? is at most ~ 0.1, and the outcomes are protoneutron stars of MjR ~ 0.1. Thus, it 
is likely that the wave amplitude is computed within ~ 10% error. The wave phase will be 
computed very accurately as indicated in [12]. 



9 



Model 


Pc {g/cw?) 




M{Mq) 


R (km) 


T/W 






(1/sCc) 


A 


A 


1 00 X 10^° 


1.503 


1.503 


2267 


8 Ql X 10~^ 


1.235 


0.994 


4.11 


oo 


B 


1.00 X 10^° 


1.485 


1.485 


1576 


5.00 X 10-=^ 


0.839 


0.994 


6.49 


0.32 


C 


1.00 X 10^° 


1.488 


1.488 


1568 


5.44 X 10-3 


0.841 


0.994 


8.45 


1/4 


D 


1.00 X 10^" 


1.500 


1.500 


1571 


1.01 X 10-2 


1.146 


0.994 


11.6 


1/4 



TABLE I. Central density, baryon rest-mass, ADM mass, equatorial circumferential radius, 
ratio of the rotational kinetic energy to the potential energy, non-dimensional angular momentum 
parameter, central value of the lapse function, angular velocity at the rotational axis, and A of 
rotating stars chosen as initial conditions for stellar core collapse simulations. 



III. INITIAL CONDITION AND COMPUTATIONAL SETTING 

Rotating stellar core in equilibrium with the F = 4/3 polytropic equation of state [see 
Eq.(lO)] is given as the initial condition for simulations. Following [10], the central density 
is chosen as pc — 10^° g/cm^ irrespective of the velocity profile. 

The velocity profiles of equilibrium rotating stellar cores are given according to a popular 
relation [24,25] 

u\ = zul{na-n), (16) 

where fla denotes the angular velocity along the rotational axis, and zud is a constant. In 
the Newtonian limit, the rotational profile is written as 

2 

^ = 2 - (17) 

Thus, Wd indicates the steepness of differential rotation. In this paper, we pick up the rigidly 
rotating case in which zud — > oo (referred to as model A) and differentially rotating cases 
with A = Wd/Re = 0.32 (referred to as model B) and 1/4 (referred to as models C and D), 
where Re is the coordinate radius at an equatorial surface. In the rigidly rotating case, we 
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chose the axial ratio of polar radius to equatorial radius as 2/3. With this choice, the angular 
velocity at the equatorial stellar surface is nearly equal to the Keplerian velocity Namely, for 
the rigidly rotating case, a rapidly rotating initial condition with nearly maximum angular 
velocity is chosen. In the differentially rotating case, we chose stars of ratio of the rotational 
kinetic energy T to the gravitational potential energy 1^ as ~ 0.005 and ~ 0.01 where 



Here, W is defined to be positive. In Table I, several quantities for the models adopted in 
the present numerical computation are summarized. 



equilibrium states with T/W ^ 0.01. With such an initial condition, the collapsing stellar 
core often forms a differentially rotating star of a highly nonspherical shape and of a high 
value of T/W [7]. It is also known that rapidly rotating neutron stars of a high degree 
of differential rotation is dynamically unstable against nonaxisymmetric deformation (e.g., 
[26] and references therein). To simulate the collapse with a high initial value olT/W, a 
nonaxisymmetric simulation will be necessary. Since our attention here is restricted to the 
axisymmetric case, we do not choose such initial conditions. 

Simulations are performed for four initial conditions hsted in Table I. Models A and B 
are almost the same initial conditions as models A1B3 and A3B2 in [10]. (Note that the 
value of A for model B is ~ 5 x 10^ km which is approximately equal to that for model 
A3B2 in [10].) Careful comparison of present numerical results with those in [10] is carried 
out using these two models. A variety of values of Fi, V2, and Fth are adopted to investigate 
the dependence of numerical results on the equations of state: Fi is chosen as 1.28, 1.30, 
1.31, and 1.32, F2 as 2 and 2.5, and Fth as 1.35, 1.5, and 5/3. The selected sets are listed in 
Table II. 

During the collapse, the central density increases from lO^'' g/cm"^ to ~ 5 x lO^'^ g/cm^. 
This implies that the characteristic length scale of the system varies by a factor of ~ 100. 




(18) 



(19) 



For the differentially rotating case with a small value of ^(< 1), it is possible to make 
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Model 


Ti 


r2 




Al 


1.32 


2.5 


1.5 


A2 


1.31 


2.5 


1.5 


A3 


1.28 


2.5 


1.5 


A4 


1.32 


2.5 


1.35 


A5 


1.32 


2.5 


5/3 


A6 


1.32 


2.0 


1.5 


Bl 


1.32 


2.5 


1.5 


B2 


1.30 


2.5 


1.5 


CI 


1.32 


2.5 


1.5 


C2 


1.30 


2.5 


1.5 


C3 


1.32 


2.0 


1.5 


C4 


1.32 


2.5 


1.35 


D 


1.32 


2.5 


1.0 



TABLE II. Selected sets of Ti, Fa, and Tth- 



One of the computational issues in a stellar core collapse simulation is to guarantee numerical 
accuracy against the significant change of the characteristic length scale. In the early phase 
of the collapse (infall phase; see Sec. IV A) in which it proceeds in a nearly homologous 
manner, we may follow the collapse with a relatively small number of grid points by moving 
the outer boundary inward while decreasing the grid spacing, without increasing the grid 
number by a large factor. As the collapse proceeds, the central region shrinks more rapidly 
than the outer region does and, hence, a better grid resolution is necessary to accurately 
follow such a rapid collapse in the central region. On the other hand, the location of the 
outer boundaries cannot be changed by a large factor to avoid discarding the matter in the 
outer envelopes. 
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Model 


Ax (initial) 


L (initial) 


Ax (final) 


L (final) 


A 


3.775 


2265 


0.4719 


1180 


A (low resolution) 


6.292 


2643 


0.7865 


1337 


B 


2.601 


1613 


0.3251 


813 


B(low resolution) 


3.902 


1639 


0.4877 


829 


C 


2.613 


1568 


0.3267 


817 


D 


2.617 


1570 


0.3271 


818 



TABLE III. The initial and final grid spacings and location of the outer boundaries along the 
X and z axes for models A-D. The unit is kilometer. 



To compute such a collapse accurately while saving the CPU time efficiently, a regrid- 
ding technique as described in [27] is adopted. The regridding is carried out whenever the 
characteristic radius of the collapsing star decreases by a factor of a few. At each regrid- 
ding, the grid spacing is decreased by a factor of 2. All the quantities in the new grid are 
calculated using the cubic interpolation. To avoid discarding the matter in the outer region, 
we also increase the grid number at the regridding, keeping a rule that the discarded baryon 
rest-mass has to be less than 3% of the total. 

Specifically, N and L in the present work are chosen in the following manner. First, we 
define a relativistic gravitational potential $c = 1 — etc (^c > 0), which is ~ 0.006 at i = for 
all the models chosen in this work. Since $c is approximately proportional to M/R, can 
be used as a measure of the characteristic length scale for the regridding. From t = to the 
time at which — 0.025, we set = 620. Note that the equatorial radius is initially covered 
by 600 grid points. At $c = 0.025, the characteristic stellar radius becomes approximately 
one fourth of the initial value. Then, the first regridding is performed; the grid spacing 
is changed to the half of the previous one and the grid number is increased to N = 1020. 
Subsequently, the value of N is chosen in the following manner; for 0.025 < $c < 0.05, we 
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set N = 1020; for 0.05 < $c < 0.1, we set N = 1700; and for 0.1 < $c, we set N = 2500 
and keep this number until the termination of the simulations since the maximum value of 
$c is at most 0.25. In this treatment, the total discarded fraction of the baryon rest- mass 
which is located outside new regridded domains is ^ 3%. 

To check the convergence of numerical results, we also perform a few simulations using 
a low grid resolution (cf. Table III). In this case, the value of N is changed as follows: for 
$c < 0.025, N = 420; for 0.025 < $c < 0.05, N = 820; for 0.05 < $c < 0.1, N = 1300; and 
for 0.1 < 1700. 

Simulations for each model with the higher grid resolution are performed for 60000-80000 
time steps. The required CPU time for one model is about 40-70 hours using 8 processors 
of FACOM VPP 5000 at the data processing center of National Astronomical Observatory 
of Japan. 

IV. NUMERICAL RESULTS 
A. Dynamics of the collapse 

1. General feature 

Figures 1-4 show the evolution of the central density (hereafter p^) and the central value 
of the lapse function (hereafter etc) for models A - D. Figures 5 and 6 are snapshots of the 
density contour curves and the velocity vectors of {v^, v^) on the y — plane for models Al 
and CI at selected time slices around which shocks are formed. 

As described in [10], rotating stellar core collapses can be divided into three phases. 
The first one is the infall phase in which the core collapse proceeds from the onset of the 
gravitational instability triggered by the sudden softening of the equation of state due to the 
reduction of the adiabatic index. During this phase, the central density (the central value of 
the lapse function) monotonically increases (decreases) until it reaches the nuclear density 
or the centrifugal force becomes strong enough to halt the collapse. The inner part of the 
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core, which collapses nearly homologously, constitutes the inner core. The duration of the 
infall phase in the present work is between about 30 and 70 msec depending mainly on the 
value of Fi (for the smaller value of Fi, the duration is shorter) as shown in [10]. We note 
that the dynamical time aX t — defined by p~^^'^ is ^ 38.7 msec. Thus, the duration may 
be written as 0.8-1.8 p'^/^. 

The second one is the bounce phase which sets in when the densities around the central 
part exceed the nuclear density pnuc, or when centrifugal forces, which become stronger 
as the collapse proceeds due to angular momentum conservation, begin to dominate over 
gravitational attraction force. At this phase, the inner core decelerates infall in about a few 
msec (~ lOp^y/^). Because of its large inertia and large kinetic energy due to the infall, 
the inner core does not settle down to a stationary state immediately but overshoots and 
bounces back, forming shocks at the outer edge of the inner core. 

The third one is the ring-down phase or the re-expansion phase. If the centrifugal force 
is sufficiently small at the time that the density of the inner core exceeds the nuclear density, 
the bounce occurs when the central density reaches ~ 2-3pnuc due to a sudden stiffening of 
the equation of state. In this case, the inner core quasiradially oscillates for about 10 msec 
and then settles down to a quasistationary state. In the outer region, on the other hand, 
shock waves propagate outward sweeping materials which infall from outer envelopes. 

If the angular momentum in the inner region is sufficiently large, the collapse is halted by 
the centrifugal force, not by the sudden stiffening of the equation of state. In this case, the 
stellar core does not settle down to a quasistationary state. Instead, it rebounds due to the 
centrifugal force and expands to be of a subnuclear density. After the maximum expansion 
is reached, the core starts collapsing again. It repeats the bounce, the expansion, and the 
collapse for many times. During each bounce, shocks are formed at the outer region of the 
core and the oscillation amplitude is damped gradually due to the shock dissipation. 

For models A-C, the centrifugal force is not strong enough to halt the collapse and, 
hence, a protoneutron star of the central density larger than the nuclear density is formed 
irrespective of the values of Fi, F2, and Fth (see Figs. 1-3). On the other hand, for model 
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D, the angular momentum is large enough to halt the collapse and to prevent the inner 
core being compact. As a result, the outcome is an oscillating star of a subnuclear density 
(see Fig. 4). Since the amplitude of the oscillation decreases gradually, it will settle to a 
rotating star of subnuclear density eventually. The adiabatic constant of this star is ^ Fi 
that is smaller than 4/3 which is the well-known critical value against gravitational collapse 
for spherical stars, and 1.329 which is an approximate critical value for rigidly rotating 
stars [28]. This indicates that the centrifugal force by a rapid and differential rotation plays 
an essential role for the stabihzation against gravitational collapse. According to [29], the 
criterion of the stability for slowly rotating stars is given by 

T M 
= 3Fi - 4 - 2 — (3Fi - 5) - fc— > 0, (20) 

where /c is a constant which is ~ 6.75 for n = 3 and T/W = [30]. For the Newtonian 
polytropes with n ~ 3, the stellar radius is given by 

/M\V3 / M \l/3/ n \-l/3 

R^2.35( — ] f«73km(^^) ( J . (21) 

\pj \1.5MqJ VlOi4g/cmV ^ ^ 

Thus, M/R will be ~ 0.03 for ~ 10^^ g/cm^. The value of T/W for dynamical stars is 
not exactly defined in general relativity, but assuming that it approximately increases as 
1/i? oc 1 — QJc for a fixed value of M, we can infer that the value of T/W would be ~ 0.15- 
0.2 for etc ~ 0.9 and k — 6.75. Therefore, Qc would be ~ 0.1-0.2, and, hence, the rotating 
star would satisfy the stability condition against gravitational collapse. On the other hand, 
the expected value of T/W is so large that the formed differentially rotating star may be 
unstable against a nonaxisymmetric deformation [26] . This suggests that to clarify the fate 
of this star, it would be necessary to perform a nonaxisymmetric simulation [18]. However, 
such a simulation is beyond scope of this paper and, hence, particular attentions are paid 
only to models A-C in this paper. 

As Figs. 1(a) and 2(a) indicate, the evolution of the central density and the central 
value of the lapse function depends strongly on the value of Fi. For the smaller value of 
Fi, the depleted pressure at i = is larger. As a result, the collapse is accelerated more 
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and the elapsed time in the infall phase is shorter. Also, since the depleted fraction of the 
pressure is larger in the central region than in the outer region, the collapse in the central 
region proceeds more rapidly. This results in a less coherent collapse for the smaller value 
of Fi. This effect makes the mass of a protoneutron star at its formation smaller and is 
reflected in the value of etc in the ring-down phase which depends on the compactness of 
the protoneutron star. On the other hand, the final value of pc depends only weakly on the 
value of Fi. This indicates that for the smaller value of Fi, the formed protoneutron star 
has a more centrally concentrated structure. 

In Figs. 1(b) and 2(b), the evolution of the central density and the central value of the 
lapse function for different values of Fth with fixed values of Fi(= 1.32) and F2(= 2.5) is 
compared. Recall that the value of Fth determines the strength of shocks at the bounce and 
at their subsequent propagation. Thus, the results here show that a moderate change of the 
value of Fth from 1.35 to 5/3 weakly modifies the evolution of the formed protoneutron stars. 
For the smaller value of Fth, the final value of the central density (central lapse) is larger 
(smaller). This is simply because the amount of matter that accretes to the protoneutron 
star increases and, hence, the compactness increases with the decrease of the value of Fth- 
For the larger value of Fth, the oscillation amplitude of Pc is larger. This is due to the fact 
that the stronger shocks result in the larger amplitude of the oscillation of the core. 

In Figs. 1(c) and 2(b), the evolution of the central density and the central value of the 
lapse function is compared for different values of F2 with fixed values of Fi(= 1.32) and 
rth(= 1-5). [Compare the solid and dottcd-dashed curves in Fig. 2(b).] Since the equation 
of state for a protoneutron star is stiffer for the larger value of F2, the maximum density 
at the bounce, the final relaxed value of pc, and the compactness of the quasistationary 
neutron star are smaller. Since the infall proceeds deeply inside the core, the amplitude of 
the oscillation for the central density in the ring-down phase is larger for the smaller value 
of F2. 

Figure 3 shows the evolution of the central density and the central value of the lapse 
function (a) for models Bl and CI and (b) for models B2 and C2. The values of Fi, F2, and 
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Fth are identical between models Bl and CI and between models B2 and C2. Furthermore, 
the values of T/W for the initial condition are approximately equal. Therefore, the difference 
of the numerical results comes from the angular velocity profile of the initial conditions. 
Figure 3 indicates that the degree of differential rotation at i = is reflected significantly in 
the oscillation and evolution of the formed protoneutron stars. The quantitative differences 
are summarized as follows: (i) the time at the bounce tb for models CI and C2 is slightly 
larger than that for models Bl and B2, respectively; (ii) the maximum value of the central 
density for models CI and C2 is slightly smaller than that for models Bl and B2, respectively; 
(iii) the amplitude of the oscillation of the central density and central value of the lapse 
function in the ring-down phase is larger for models CI and C2. The results (i) and (ii) are 
simply due to the fact that the centrifugal force around the central region for models CI 
and C2 is slightly larger and plays a stronger role for halting the collapse. The result (iii) 
indicates that the small increase of the angular velocity around the central region in the 
initial condition can significantly modify the evolution of the central density. All the results 
(i)-(iii) also show that the oscillation of the central density of the formed protoneutron stars 
depends strongly on the initial angular velocity profile. 

Effects of differential rotation of the initial condition are also refiected significantly in the 
shape of the formed protoneutron stars. In the collapse of a rigidly rotating progenitor, the 
formed protoneutron star has a slightly nonspherical shape (see Fig. 5). On the other hand, 
in the collapse of a differentially rotating progenitor, a protoneutron star of a fiattened and 
nonspherical shape is the outcome (see Fig. 6). This difference results from the fact that 
the inner region is more rapidly rotating in the case of the differentially rotating progenitor. 
It is worthy to note that the value of T/W for model A is about 1.6 times as large as that 
for model C. However, the angular velocity at the rotational axis for model A is about half 
of that for model C. Thus, T/W alone is not a good indicator for measuring the significance 
of the centrifugal force in rotating stellar core collapses (nor is the nondimensional angular 
momentum parameter J/M^). Obviously, the local distribution of the angular momentum 
plays a more important role for determining the shapes of the formed protoneutron star and 
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shocks. 

Convergence of the numerical results is achieved well in the present computation. In 
Figs. 1(d), 3(a), and 3(b), we show the numerical results with a low grid resolution for 
models Al, Bl, and B2 (dotted curves). It is found that the evolution of the central density 
and the central lapse in the low-resolution simulation agrees with that in the high-resolution 
one within a small error (except for the very late time for which the numerical error seems to 
be accumulated for the low- resolution simulation). This indicates that the grid resolutions 
adopted in the present numerical simulation are fine enough to yield a convergent numerical 
result. 

2. Comparison with a previous work 

Here, we compare the numerical results for models Al, A2, A3, Bl and B2 with those 
for models A1B3G2, A1B3G3, A1B3G5, A3B2G2, and A3B2G4 in [10], respectively For 
these models, both groups adopt the almost identical initial conditions. 

Table IV shows the time at a maximum density achieved, the maximum density, and 
the maximum amplitude of gravitational waves for the numerical results computed by two 
groups. In Fig. 7, we also compare the evolution of the central density. It is found that the 
numerical results by two groups agree within a small error both for models A and B. Only for 
model Bl, the time at the maximum density achieved slightly disagrees with that for A3B2G2 
by ~ 0.3 msec, but besides this disagreement, the shape of pc as a function of time agrees 
well each other even in this case. Recall that in [10], the conformal flatness approximation 
to the Einstein equation is adopted while our results are fully general relativistic. This 
indicates that the conformal flatness approximation is a good approximate formulation of 
general relativity for computing axisymmetric rotating stellar core collapses to a neutron 
star. 

In a precise comparison, the following small systematic disagreements between two results 
should be also addressed: (i) the maximum density achieved in our results is slightly larger 
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for model A and slightly smaller for model B; (ii) the time at the maximum density is slightly 
delayed in our results, and this tendency is stronger for the larger value of Fi (i.e., for the 
longer infall time); (iii) for the larger value of Fi, the central density in the relaxed final 
stage is slightly smaller in our results. 

It is difficult to specify the particular reason for these disagreements. There are several 
plausible candidates. First, computational settings are different between two groups. In our 
simulation, we adopted a uniform grid changing the grid spacing and grid number, while 
in [10], 200 radial grid points with the logarithmic grid spacing were taken throughout the 
simulation. In our case, the grid spacing is smaller than 0.5 km in the bounce and ring- 
down phases, although it is larger than 0.5 km in the infall phase. On the other hand, the 
minimum grid spacing is about 0.5 km in [10] for all the phases. These differences may 
yield the disagreements. Actually, we find that varying the grid resolution results in a small 
change of the time at the maximum density achieved for models Al and Bl [cf. Figs. 1(d) 
and 3(a)]. Secondly, the slicing condition is slightly different between two groups. In [10], 
the maximal slicing condition, K^!' = 0, was adopted, while in our numerical simulation, 
the condition is only approximately satisfied [16]: The equations K^^ = = dtK^^ lead 
to an elliptic-type equation for a. In the exact maximal slicing condition, this equation is 
iteratively solved until a convergence is achieved. In our case, we stop the iteration before 
the complete convergence is achieved to save the computational time. Thus, ^ 0. This 
difference may result in a systematic deviation of the coordinate time at the maximum 
density. Thirdly, the initial conditions adopted by two groups are not completely identical, 
since the equilibrium rotating stars for the initial conditions are computed with different 
numerical implementations. The values of T/W and A may well have disagreement of 
magnitude ^ 1%. This may affect the subsequent numerical evolution slightly. 

On the other hand, the difference of the adopted formulations for the gravitational field 
is unlikely to be the reason for the disagreement. This is because the deviation of the 
conformal metric from Sij is very small (typical absolute magnitude is of order ~ 10"^ 
for each component) in our numerical results. Therefore, we infer that the magnitude of 
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Model 


tb 


Pmax (g/cm'^) 


{rh_^_ )max (cm) 


Al 


69.5 


4.12 


561 


A1B3G2 


69.5 


4.02 


469 


A2 


48.7 


4.28 


215 


A1B3G3 


48.6 


4.23 


180 


A3 


30.3 


4.98 


32.7 


A1B3G5 


30.2 


4.55 


33.9 


Bl 


69 8 


3 93 


731 


A3B2G2 


69.5 


4.10 


596 


B2 


39.3 


3.92 


182 


A3B2G4 


39.3 


4.05 


141 



TABLE IV. Comparison between the present (upper) and previous numerical results by Dim- 

mclmcicr et al. (lower). The time at bounce, the maximum density achieved, and the maximum 
amplitude of gravitational waves are shown for two numerical results. 



the systematic error due to the conformally flatness approximation seems to be smaller than 
that due to other reasons. 

B. Gravitational waveforms 

1. General feature 

Gravitational waveforms are computed in terms of the quadrupole formula described in 
Sec. II C. Since fully general relativistic simulations are performed, gravitational waves 
should be computed from the metric in a wave zone. However, we have found that it is 
not possible, since the amplitude is smaller than the numerical noise. An estimate by the 
quadrupole formula indicates that the maximum amplitude of gravitational waves is smaller 
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than 10~^ in the local wave zone for r ~ A where A denotes the wave length which is typically 
several hundred km. 

As illustrated in a previous paper [12], approximate gravitational waveforms can be com- 
puted in terms of a quadrupole formula for highly relativistic, highly oscillating, and rapidly 
rotating neutron stars. In rotating stellar core collapses to a neutron star, gravitational 
waves are dominantly emitted during the bounce and ring-down phases. Such gravitational 
waves are excited by oscillations of a formed protoneutron star. Thus, it is likely that the 
present approach can yield high-quality approximate gravitational waveforms besides possi- 
ble underestimation of the amplitude by ~ 10% due to absence of higher general relativistic 
corrections. 

Figures 8(a)-(d) show gravitational waveforms for model A with various sets of Fi, F2, 
and Fth- The waveforms for models Al, A2, A4, A5, and A6 are classified into type I 
according to Dimmelmeier et al. [10]. Properties of the type I gravitational waveforms 
can be summarized as follows: During the infall phase, a precursor whose amplitude and 
characteristic frequency increase monotonically with time is emitted due to the infall and 
the flattening of the rotating core. The duration of the infall phase is ^ 40 msec and longer 
than a dynamical time scale defined at i = as p~^/^ ~ 40 msec. In the bounce phase, 
spiky burst waves are emitted for a short time scale ~ 1 msec, and the amplitude and the 
frequency of gravitational waves become maximum. In the ring-down phase, gravitational 
waves associated with several oscillation modes of a formed protoneutron star are emitted 
and its amplitude is gradually damped due to shock dissipation at the outer edge of the 
protoneutron star. 

For model A3 [cf. Fig. 8(c)] for which the simulation is performed with a small value 
of Fi as 1.28, the waveforms are qualitatively different from those for others: A sharp and 
distinguishable peak is not found at the bounce. Soon after the precursor emitted during 
the infall phase, the ring-down waveforms appear to be excited. An outstanding feature is 
that the amplitude in this case is much smaller than that for Fi = 1.31 and 1.32 although 
the wave length is not significantly different from those for other models. According to [10], 
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this type of the waveforms is classified into type III. 

In Fig. 8(a), the waveforms for models Al, A4, and A5 are presented. For these models, 
we adopt Fi = 1.32 and F2 = 2.5, so that only the value of Fth is different. In the infall 
phase, the waveforms for three models are very similar. This is natural because as long 
as the density is smaller than pnuc; the magnitude of Pth is much smaller than that of the 
cold part. Clear differences in the wave phase, wave length, and amplitude are observed 
in the bounce and ring-down phases. The reasons for them are explained as follows: The 
smaller magnitude of Pth results in the slightly shorter infall time as reflected in the time at 
which the amplitude becomes maximum. As a consequence, the difference of the wave phase 
is yielded. Stronger shock heating, which generates larger thermal energy, also results in 
smaller compactness of the formed protoneutron stars. This leads to the results that for the 
larger value of Fth, the gravitational wave length, which in general increases with the stellar 
radius for a given mass, becomes longer, and the amplitude, which is larger for stronger 
shock heating, is larger. 

Slight change of the value of Fi, which determines the dynamics of the infall phase, 
signiflcantly modifles gravitational waveforms. Comparison among Figs. 8(a)-(c) clarifles 
that with the decrease of the value of Fi, the amplitude of gravitational waves decreases 
systematically. The reason for this is explained as follows: For the smaller value of Fi, the 
central region collapses more rapidly than the outer region does. This results in a smaller 
core mass at the bounce for the smaller value of Fi. The amphtude of gravitational waves 
increases with the increase of the core mass for a fixed value of the density and, therefore, 
it is smaller for the smaller value of Fi. 

In Fig. 8(d), we compare the waveforms of different values of F2 with fixed values of Fi 
and Fth. It is found that the difference of the waveforms between two models appears only in 
the bounce and ring-down phases. This is natural because the value of F2 does not affect the 
infall phase and mainly determines the equations of state and the radius (or compactness) 
of the formed protoneutron stars. Recall that the smaller value of F2 results in the larger 
compactness of the protoneutron star. This fact is reflected in slightly shorter wave length 
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and larger amplitude of gravitational waves in the ring-down phase for the smaller value of 

r2. 

Figure 9 displays gravitational waveforms for model C. As in the case of model A, the 
waveforms are divided into three parts (precursor, spike, and ring-down), but the qualitative 
feature of the ring-down waveforms between models A and C is different. For example, 
compare the waveforms for models Al and CI for which the values of Fi, F2, and Fth are 
identical. For model Al, the waveforms are modulated only in the early ring-down phase 
(e.g., for t ~ 70-73 msec). In the late ring-down phase (e.g., for i > 73 msec for model Al), 
they are fairly periodic and appear to be composed mainly of one or two eigen oscillation 
modes of the formed protoneutron star. On the other hand, for model CI, the waveforms are 
not very periodic and highly modulated throughout the ring-down phase. In this case, several 
eigen modes of a formed protoneutron star appear to constitute gravitational waveforms. 
Such modulated waveforms are likely to be due to the fact that the formed protoneutron star 
is rapidly and differentially rotating and the oscillation modes are excited in a complicated 
manner at the bounce. 

In Fig. 9(a), we compare the waveforms of different values of Fth with fixed values of Fi 
and F2. As in Fig. 8(a), for the smaller value of Fth, the maximum amplitude is reached 
at an earlier time, the wave length during the bounce and ring-down phases is longer, and 
the amplitude is smaller. These are universal features independent of the initial rotational 
velocity profiles. However, in contrast to Fig. 8(a), the waveforms in the ring-down phase 
for models CI and C4 are not very similar. Thus, small change of Fth from 1.35 to 1.5 
significantly modifies the ring-down waveform in the case of differentially rotating initial 
velocity profiles. 

In Fig. 9(b), we compare the waveforms of different values of F2 with fixed values of Fi 
and Fth. In contrast to Fig. 8(d), the maximum amplitude of gravitational waves is nearly 
identical for two models. This suggests that in halting the infall, the centrifugal forces may 
play an important role to hide the effects of the difference in the value of F2. The difference 
of the ring-down waveforms between two models is qualitatively the same as that found in 
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Fig. 8(d): For the smaller value of the wave length and amphtude of gravitational waves 
in the ring-down phase are slightly shorter and larger, respectively. 

In Fig. 9(c), the waveform for model C2 is displayed. This should be compared with 
the sohd curve in Fig. 9(a) [or 9(b)] for model CI of a different value of Fi. Comparison 
between two waveforms shows that with the decrease of the value of Fi, the wave amplitude 
at the bounce and ring-down phases decreases. This property agrees with that found for 
model A and is likely to be independent of the initial rotational velocity profiles. 

To see the effect of the slight change of differential rotation parameter A, we compare 
the waveforms of models Bl (solid curve) and CI (dotted-dashed curve) in Fig. 9(d). Two 
waveforms are qualitatively similar, but for model CI the amplitude is larger and the mod- 
ulation of the amplitude is induced more. This illustrates that with a slight modification 
of the initial rotational velocity profile, the resulting gravitational waveforms are modified 
significantly. 

Figure 10 shows the gravitational waveform for model D. In this model, the collapse 
does not lead to a quasistationary protoneutron star of Pc > Pnuc- Instead, a quasiradially 
oscillating star of subnuclear density is formed and, therefore, quasiperiodic waves of a long 
period ~ 10 msec are emitted. According to [10], this is classified into the type II waveform. 

Convergence of the numerical results appears to be achieved. In Fig. 11, we display 
the numerical results with high and low grid resolutions for models Al and Bl. The grid 
spacing in the low grid resolution is about 5/3 as large as that in the high case. It is found 
that the computed gravitational waveforms depend only weakly on the grid resolution in 
our choice of the grid spacing. We conclude that the grid resolution we choose in this work 
is fine enough to compute convergent gravitational waveforms. 

2. Comparison with a previous work 

Here, we compare the gravitational waveforms computed in this paper with those in [10] 
for models Al, A2, A3, Bl, and B2. Figures 12 and 13 show the gravitational waveforms 
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computed by us (solid curves) and by Dimmelmeier et al. [10] (dashed curves). It is found 
that the waveforms in the infall phase agree very well each other. In the bounce phase, on 
the other hand, the amplitude of our results is larger than that in [10] by ~ 20% for models 
Al, A2, Bl, and B2, although they still agree qualitatively. The disagreement is outstanding 
in the ring-down phase. The amplitudes of gravitational waves in the ring-down phase for 
models Al, A2, Bl, and B2 are larger than those in [10] by a factor of ~ 2. Moreover, in 
our results, the oscillations with a nearly constant amplitude continue for several oscillation 
periods (^ 10 msec). This is not the case in the results of [10] in which the amplitude is 
damped within several msec. 

This would be partly due to the difference in the grid resolution or in the slicing condition 
adopted by two groups as mentioned in the previous section. However, the main reason is 
likely that quadrupole formulas adopted by two groups are not identical. In the quadrupole 
formula we adopt, a quadrupole moment is simply defined using a weighted rest- mass density 
and then the second time derivative is taken with no approximation. In [10], on the other 
hand, the quadrupole moment is defined using p and, in addition, when taking the second 
time derivative, they discard higher relativistic terms keeping only the lowest order post- 
Newtonian terms. 

This disagreement raises a question: what is a good quadrupole formula in general 
relativistic simulations ? An excellent quadrupole formula should yield a high-quality ap- 
proximate waveform for the true one computed from the metric in the wave zone. Thus, 
to answer the question, it is necessary to compare the gravitational waveforms computed 
by a quadrupole formula with those extracted from the metric. In [12], we calibrated the 
waveforms by performing simulations for highly relativistic, highly oscillating, and rapidly 
rotating neutron stars with M/ R ~ 0.2 and v/c ^ 0.3 and found that our quadrupole formula 
yields well- approximated waveforms; the wave phases agree well with those computed from 
the metric and the wave amplitude is computed within an error of magnitude of 0{M/R) 
or 0{v'^/c^). We believe that the waveforms presented in this paper are well-approximated 
ones in phase and within ~ 10% error in amplitude. On the other hand, the quadrupole 
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formula adopted in [10] has not been calibrated, since they adopted the conformal flatness 
approximation in which gravitational waves cannot be extracted from the metric. Thus, 
it is not clear how good their quadrupole formula is. Since the amplitudes computed by 
our quadrupole formula are underestimated by ~ 10% and the amplitudes computed in [10] 
are smaller than ours, gravitational waveforms presented in [10] may contain an error of 
magnitude more than 10-20%. 

V. SUMMARY 

We performed axisymmetric numerical simulations of rotating stellar core collapses to a 
neutron star in full general relativity, paying particular attention to gravitational waveforms 
and to comparison of our results with previous results [10]. The Einstein field equations 
are solved in the Cartesian coordinates imposing an axisymmetric condition by the Cartoon 
method [22]. The hydrodynamic equations are solved in the cylindrical coordinates (with the 
Cartesian coordinates restricting on the y — plane) using a high- resolution shock-capturing 
scheme with the maximum grid size (2500, 2500). A parametric equation of state is adopted 
to model collapsing stellar cores and formed protoneutron stars following Dimmelmeier et 
al. [10]. Gravitational waveforms are computed using a quadrupole formula proposed in [12]. 

We choose moderately rapidly rotating stars as the initial conditions for which the value 
ofT/W is between 0.005 and 0.01. Simulations are performed changing three parameters (Fi, 

and Fth) which characterize the equation of state. The dynamics of the collapse depends 
on the three parameters as well asT/W and A of the initial condition. The dependence of 
the evolution of the system and gravitational waveforms on these five parameters is studied. 
The value of Fi mainly determines the duration of the infall phase and the coherence of the 
early phase of the collapse. For the smaller value of Fi, the infall time becomes shorter and 
the collapse is accelerated more in the central region. This results in that the core mass at 
the bounce is smaller and that the magnitude of $c (which may be regarded as the depth 
of the gravitational potential) at the bounce is smaller for the smaller value of Fi. The 
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amplitude of gravitational waves becomes also smaller for the smaller value of Fi. 

The value of r2 determines the equation of state for formed protoneutron stars. Thus, 
it does not affect the evolution during the infall phase. It determines the final value of the 
central density of the formed protoneutron star and the gravitational waveforms emitted 
during the ring-down phase in which the eigen oscillation modes of the protoneutron stars 
are excited. The value of Fth determines the strength of shock waves. We choose this value 
as 1.35, 1.5, and 5/3 extending the work by Dimmelmeier et al. [10]. It is found that for 
the smaller value of this parameter, the shock heating becomes weaker and the amplitude 
of gravitational waves smaller. 

The values oiT/W and A play a significant role in determining the dynamics of collapse 
and the corresponding gravitational waveforms in particular in the bounce and ring-down 
phases. For the rigidly rotating case (^4 — oo), the maximum value of T/W is ~ 0.009 
which we choose in this paper. Even in this maximum case, the collapse leads to a neutron 
star irrespective of the values of Fi, F2, and Fth. This indicates that for rigidly rotating 
initial conditions, the neutron stars are formed soon after the collapse, irrespective of the 
angular velocity of the initial condition, with our choice of the equations of state. For the 
differentially rotating case with A — 1/4, the collapse does not lead to a neutron star but an 
oscillating star of subnuclear density is formed for T/W ^ 0.01 since the centrifugal force 
is strong enough near the rotational axis. As shown in [10], more rapidly rotating initial 
conditions with T/W » 0.01 may be constructed. For such high values of T/W, a neutron 
star will not be formed soon after the collapse. 

With a slight change of A from 0.25 to 0.32 for the initial condition, the angular velocity at 
the rotational axis is changed by a large factor even if T/W is approximately identical. As a 
result of this change, the subsequent evolution of the collapse and gravitational waveforms in 
the bounce and ring-down phases are modified significantly. This implies that the dynamics 
rotating stellar core collapses and the corresponding gravitational waveforms are sensitive 
not only to the equation of state but also to the initial angular velocity profile. 

Several simulations are performed setting the same initial conditions as those adopted in 
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[10]. It is found that the dynamics of the collapse and the bounce for such initial conditions 
are very similar to those found in [10], in which an approximate general relativistic gravity 
(the conformal flatness approximation) is assumed. This indicates that such approximate 
relativistic formulation is appropriate for computing axisymmetric rotating stellar core col- 
lapses and the subsequent formation of protoneutron stars. (Note that this is the conclusion 
for the formation of neutron stars. This may not be the case for the black hole formation.) 

Gravitational waveforms are compared with a previous result [10]. It is found that 
the waveforms are qualitatively in good agreement but not quantitatively with those in 
[10]. Either of two plausible elements could explain this disagreement. One is that the 
grid resolution and computational setting are different between two groups. This could 
modify the waveforms slightly. However, the main reason seems to be that the adopted 
quadrupole formulas by two groups are different. As mentioned in the previous section, 
there is no unique definition of the quadrupole formula for dynamical spacetimes in general 
relativity. This implies that when one attempts to use a quadrupole formula in a relativistic 
simulation, one needs to calibrate the formula in advance by performing a fully general 
relativistic simulation and by comparing the waveforms computed by the quadrupole formula 
with those computed from the metric. The quadrupole formula adopted in our study has 
been calibrated in simulations for highly relativistic, highly oscillating, and rapidly rotating 
neutron stars [12]. Thus, we believe that the quadrupole formula adopted in this paper is 
appropriate and that the numerical results presented here are approximate solutions of high 
quality. 

In this paper, we have focused on the neutron star formation and on the comparison with 
a previous work [10]. If more massive progenitor is chosen as the initial condition, a black 
hole instead of a neutron star will be formed. Formation of black holes and corresponding 
gravitational waves have been studied by several groups [31,32,17]. However, the initial 
conditions in their previous works are not very realistic for modeling a rotating stellar core 
collapse in nature. Namely, the stellar core collapse to a black hole from a realistic initial 
condition in fully general relativistic simulation is an unsolved issue. We are currently 
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working in this subject and will present the numerical results in a subsequent paper. 
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FIG. 1. Evolution of the central density and the central value of the lapse function for model 

A. (a) Models Al (solid curve), A2 (dashed curve), and A3 (dotted-dashed curve); (b) Models Al 

(solid curve), A4 (dotted curve), and A5 (dotted-dashed curve); (c) Models Al (solid curve) and 

A6 (dotted curve); (d) Model Al with high (solid curve) and low grid resolutions (dotted curve). 
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FIG. 2. (a) The same as Fig. 1 but for models CI (solid curves), C2 (dashed curves), and C3 

(dotted-dashed curves), (b) The same as Fig. 1 but for models CI (solid curves), C3 (dotted-dashed 

curves), and C4 (dotted curves). 
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FIG. 3. (a) The same as Fig. 1 but for model Bl with high (solid curve) and low grid resolutions 

(dotted curve). For comparison, the results for model CI (dashed curve) are also shown, (b) The 

same as (a) but for model B2 with high (solid curve) and low grid resolutions (dotted curve) and 

for model C2 (dashed curve). 
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FIG. 4. Evolution of the central density and the central value of the lapse function for model D. 
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FIG. 5. Snapshots of the density contour curves of p and of the velocity field of (v^ ,v^) at 
selected time slices around which shocks are formed for model Al. The contour curves are drawn 
for p/pnuc = 3 X 10-°-^^ with j = 0, 1, 2, • • • , 15. 
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FIG. 7. Comparison between the evolution of the central density computed in this paper (solid 

curves) and in Dimmelmeier et al. (dashed curves) (a) for models A1-A3 and (b) for models Bl 

and B2. 
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FIG. 8. Gravitational waveforms for model A; (a) Models Al (solid curve), A4 (dotted-dashed 

curve), and A5 (dotted curve); (b) Model A2; (c) Model A3; (d) Models A6 (solid curve) and Al 
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FIG. 9. Gravitational waveforms for model C; (a) Models CI (solid curves) and C4 (dashed 

curves); (b) Models CI (solid curves) and C3 (long-dashed curves); (c) Model C2; (d) Models Bl 

(solid curve) and CI (dotted-dashed curve). 
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FIG. 10. Gravitational waveforms for model D. 
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FIG. 11. Gravitational waveforms (a) for model Al and (b) for model Bl with high (solid 
curve) and low grid resolutions (dashed curve). 
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FIG. 12. Comparison between gravitational waveforms computed in this paper (solid curves) 
and in Dimmelmeier et al. (dashed curves) for models A1-A3 [(a)-(c)]. 
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